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Hybrid  DC  power  sources  which  consist  of  fuel  cells,  photovoltaic  and  lithium-ion  batteries  provide  clean, 
high  efficiency  power  supply.  This  hybrid  DC  power  sources  can  be  used  in  many  applications.  In  this 
work,  a  model-based  fault  detection  methodology  for  this  hybrid  DC  power  sources  is  presented.  Firstly, 
the  dynamic  models  of  fuel  cells,  photovoltaic  and  lithium-ion  batteries  are  built.  The  state  space  model 
of  hybrid  DC  power  sources  is  obtained  by  linearizing  these  dynamic  models  in  operation  points.  Based 
on  this  state  space  model  the  fault  detection  methodology  is  proposed.  Simulation  results  show  that 
model-based  fault  detection  methodology  can  find  the  fault  on  line,  improve  the  generation  time  and 
avoid  permanent  damage  to  the  equipment. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Currently,  most  of  the  energy  demand  in  the  world  is  met  by 
fossil  and  nuclear  power  plants.  A  small  part  is  drawn  from  renew¬ 
able  energy  technologies  such  as  wind,  solar,  fuel  cell,  biomass  and 
geothermal  energy.  Solar  energy,  fuel  cells  and  wind  energy  have 
experienced  a  remarkably  rapid  growth  in  the  past  1 0  years  because 
they  are  pollution-free  power  sources.  Additionally,  they  generate 
power  near  the  load  centers,  which  eliminates  the  need  to  run  high- 
voltage  transmission  lines  through  rural  and  urban  landscapes. 
Hybrid  DC  power  sources  which  consist  of  fuel  cells  (FC),  photo¬ 
voltaic  (PV)  and  lithium-ion  batteries  provide  clean,  high  efficiency 
power  supply.  This  hybrid  DC  power  sources  can  be  used  in  many 
applications  such  as  home  power  supply,  electrical  car  and  satel¬ 
lite  power  systems.  But  the  structure  of  hybrid  DC  power  sources 
is  complex  and  vulnerable  to  faults.  For  examples,  fuel  cells  have 
faults  in  the  air-reaction  blower,  the  refrigeration  system,  hydrogen 
pressure  and  increase  of  fuel  crossover  [1,2].  The  energy  losses  of 
PV  system  may  be  resulted  from  four  different  fault  categories:  sus¬ 
tained  zero  efficiency  faults;  brief  zero  efficiency  faults;  shading; 
and  non-zero  efficiency  non-shading  faults  [3].  So  it  is  necessary  to 
study  fault  detection  and  diagnosis  of  this  hybrid  DC  power  sources. 

Fault  detection  is  a  subfield  of  control  engineering  which 
concerns  itself  with  monitoring  a  system,  identifying  when  a 
fault  has  occurred  and  pinpoint  the  type  of  fault  and  its  loca¬ 
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tion.  A  fault  is  defined  as  an  unexpected  change  in  a  system 
with  component  malfunction  or  variation  in  operating  condi¬ 
tion.  Faults  in  a  dynamic  system  can  take  many  forms,  such  as 
actuator  faults,  sensor  faults  and  abrupt  changes  of  some  param¬ 
eters.  The  representative  fault  detection  approaches  are  model 
free  method,  knowledge-based  method  and  model-based  method. 
The  model-based  fault  detection  approaches  include  parity  space- 
based  approach  [4],  eigenstructure  assignment-based  approach  [5], 
parameter  identification-based  approach  [6]  and  observer-based 
approach.  Among  these  methods,  observer-based  fault  detection  is 
one  of  the  most  effective  methods  and  has  obtained  much  more 
attention. 

In  this  paper,  a  model  based  fault  detection  method  is  proposed 
to  execute  on-line  fault  diagnosis  of  hybrid  DC  power  sources.  The 
structure  of  this  paper  is  the  following:  in  Section  2,  the  models  of 
the  hybrid  DC  power  sources  were  built.  In  Section  3,  the  proposed 
model-based  fault  diagnosis  methodology  was  described.  In  Sec¬ 
tion  4,  the  hybrid  DC  power  sources  system  is  used  to  illustrate  the 
proposed  fault  diagnosis  methodology  with  the  fault  scenarios  that 
can  appear.  Finally,  in  Section  5,  conclusions  are  drawn,  based  on 
the  results  from  the  numerical  analyses  in  Section  4. 

2.  Model  of  hybrid  DC  power  sources 

2. 1 .  The  structure  of  hybrid  DC  power  sources 

The  hybrid  DC  power  system  consists  of  a  photovoltaic  system, 
a  PEM  fuel  cell  power  source,  and  a  lithium-ion  battery,  which  are 
connected  to  the  same  DC  voltage  bus  through  appropriate  DC/DC 
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DC/DC 


Fig.  1.  Structure  of  hybrid  power  sources. 


power  converters  and  controller.  Fig.  1  illustrates  the  structure  of 
the  proposed  hybrid  DC  power  system.  There  are  two  main  power 
sources:  PV  and  FC.  When  the  supplies  of  FC  and  PV  excess  the  load 
demands  the  lithium-ion  battery  is  a  power  source.  In  the  opera¬ 
tion  of  hybrid  DC  power  sources  the  PV  provides  as  much  power  as 
possible  to  the  load.  In  order  to  operate  in  maximum  power  point,  a 
maximum  power  point  tracker  (MPPT)  controller  must  be  installed 
in  PV  system.  The  voltage  of  DC  bus  is  288  V  and  the  voltage  of  PV 
system  is  varied  between  10  V  and  100  V,  it  is  necessary  to  insert  a 
DC/DC  converter  between  the  PV  system  and  DC  bus.  The  FC  is  to 
supply  the  rest  power  of  load.  A  DC/DC  converter  is  also  installed 
between  the  FC  and  DC  bus.  The  lithium-ion  battery  supplies  tran¬ 
sient  power  to  peak  load  demands  or  absorbs  transient  power  from 
the  main  sources.  The  currents  and  voltages  of  the  PV  panel,  the  fuel 
cell,  and  the  lithium-ion  battery  are  monitored.  The  measured  volt¬ 
ages  and  currents  are  fed  into  the  power  management  subsystem, 
which  is  used  to  controls  the  power  flowing  from  each  source  of 
energy,  and  allocates  the  available  power  to  recharge  the  battery  if 
possible. 

2.2.  Model  of  PEM  fuel  cell 

We  build  water  balance  equation  by  investigating  stack  voltage 
model,  water  distribution  of  membrane  and  water  transport  pro¬ 
cess  in  cathode  and  anode.  In  this  paper  direct  injection  of  water 
into  cathode  is  used  for  humidification.  The  following  equations 
are  stack  voltage  model,  water  distribution  of  membrane  and  water 
balance  equation  of  cathode  and  anode. 

We  assume  that  single  fuel  cell  is  symmetrical  the  stack  voltage 
can  be  represented  by  multiplying  the  voltage  of  single  fuel  cell 
with  number  of  fuel  cell. 

Vst  =  n  ■  Vfc  (1) 

The  voltage  of  single  fuel  cell  VfC  is  calculated  as  a  function  of 
current,  cathode  pressure,  reactant  partial  pressures,  stack  tem¬ 
perature  and  membrane  humidity.  There  are  a  lot  of  voltage  model 
in  the  literature,  we  adopt  the  semi-empirical  model  in  this  paper 
[7]: 

Vfc  =  Eq  —  V act  —  V 0\ im  —  Vconc  (2) 

Eq  is  the  thermodynamic  potential  of  the  cell  representing  its 
reversible  voltage: 

E0  =  1.229  -  0.85  x  10“3(T  -  298.15) 

+  4.3085  x  10“5r  [ln(pH2)  +  0.5  ln(po2)J  (3) 

where  T is  operating  temperature,  Ph2  is  the  hydrogen  pressure  and 
Po2  is  the  oxygen  pressure. 


Activation  voltage  loss  Vact  is  the  voltage  drop  due  to  the  acti¬ 
vation  of  the  anode  and  the  cathode: 

Vact  =  V0  +  Va(  l-e-Cl1')  (4) 

where  C\  is  parameter  [7,8],  Va  and  V0  can  be  calculated  as  follows: 
V0  =  0.279-8.5  x  10-4(T-  298.15)  +  4.308  x  lO-5! 


In 


( Pea  ~  Psat  \  1  |  ( 

lT0i325  j +  2lnl 


0.01 173(pca  —  Psat) 
1.01325 


)] 


(5) 


Va  =  (-1.618  X  10“5T  + 1.618  X  10"2)(q^ 
+  (1.8  x  10“4T-0.166)  (2^ 


173 


+  Psat 


1173 


+  Psat 


+  (-5.8  x  10“4r  + 0.5736)  (6) 

Vohm  is  the  ohmic  voltage  drop  associated  with  the  conduction  of 
protons  through  the  electrolyte,  and  of  electrons  through  the  inter¬ 
nal  electronic  resistance  and  calculated  as  follows: 

trn 


Vohm  ~  i  '  Rohm  —  *  ' 


(7m 


(7) 


where  tm  is  the  thickness  of  membrane,  om  is  the  membrane  con¬ 
ductivity  and  calculated  as  follows: 


om=b,  exp  (ft,  (2^-1)) 


(8) 


in  which  b2  constant  [9],  hi  function  of  water  content  of  membrane 

Am. 

Vconc  is  the  concentration  voltage  loss  resulted  from  the  drop  in 
concentration  of  the  reactants  as  they  are  consumed  in  the  reaction: 

Vconc  =  i  (  t——  )  (9) 

where  c2,  c3  and  imax  are  constants  [8]  and  depend  on  stack  tem¬ 
perature  and  pressure  of  reactants. 

The  governing  equation  for  hydrogen  and  water  in  the  anode 
can  be  written  as: 
drriHn 

-ff-  =  wti2jn  -  wH2,0Ut  -  wH2,rec 

dm  w,  an 


(10) 


^  ~  ^w,an,z'n  WWjan,out  Ww,mbr 


(11) 


The  hydrogen  and  water  flowing  to  the  anode  are  calculated  by 

WH2,in  =  (l-a)Wan,in  (12) 

V/W,an,in  —  ^^4fn,z'n  (18) 

where  a  is  inlet  flow  humidity  ratio  in  anode. 

The  outlet  hydrogen  and  vapor  mass  flow  rate  are  calculated  by 
the  following  equations: 

V/an,out 


VV\\  2, an, out  — 
Ww,an,out  — 
(O  an,  out  — 


1  +  (Oan,out 
( Oan,out 


W, 


1  +  0)an,out 

Mv  m\\2R\\2 


an, out 


(14) 

(15) 

(16) 


Mh2  mv,anRH2 

Inside  the  cathode  there  are  three  elements:  oxygen,  nitrogen 
and  water.  The  balance  equation  for  the  mass  of  these  elements  can 
be  written  as: 
dm02 


~  W02 , in  V/q 2,out  V/q2  rec 

-^  =  wN2jn-wN2,0Ut 


(17) 

(18) 
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— —  —  WW  Ca,in  ~  Ww,Ca,out  +  W w,gen  +  Ww  mbr  +  ^inj  (19) 

where  Wo2jin  is  the  mass  flow  rate  of  oxygen  gas  entering  the  cath¬ 
ode  ;  W02 1  out  is  the  mass  flow  rate  of  oxygen  gas  leaving  the  cathode ; 
W0 2jec  is  the  mass  flow  rate  of  oxygen  reacted;  WNl  in  is  the  mass 
flow  rate  of  nitrogen  gas  entering  the  cathode;  W^2t0Ut  is  the  mass 
flow  rate  of  nitrogen  gas  leaving  the  cathode;  Ww  ca  in  is  the  mass 
flow  rate  of  vapor  entering  the  cathode;  WW!ca,out  is  the  mass  flow 
rate  of  vapor  leaving  the  cathode;  WW5gen  is  the  mass  flow  rate  of 
vapor  generated  in  the  fuel  cell  reaction;  Ww  inj  is  the  mass  flow 
rate  of  injected  water  from  humidifier;  Ww,mem  is  the  mass  flow 
rate  of  water  transfer  across  the  membrane. 

Inlet  mass  of  oxygen,  nitrogen  and  vapor  is  related  to  mass  flow 
rate  and  humidity  of  inlet  air: 


w  1 

02,m— yo  2(l+^atm)Wcajn 

(20) 

VV",— yN  2(1  +  i2atm)Wc0.n 

(21) 

w  ■  -  Q atm 

W’ca’m  (1  +C2atmWca,in 

(22) 

where  £2atm  is  the  humidity  ratio  of  air,  y02  and  yN2 
mole  fraction  and  nitrogen  mole  fraction  of  dry  air. 

are  oxygen 

mo2 

y°2  ~  X°2  jyjatm 

(23) 

v  n  y  i  Mn2 

Yn2  (1  x02)Matm 

(24) 

Xq2  is  volume  content  of  oxygen  in  dry  air. 

M%tm  =  Xq2Mq2  +  (1  -  Xq2  )Mn2 

(25) 

The  calculations  ofW02,out,  WN2,0Ut,  Ww,ca,out,  Ww,mem  can  find 
in  Ref.  [8]. 

The  rate  of  oxygen  consumed  in  the  reaction  and  water  gener¬ 
ated  in  the  reaction  is  calculated  by 

M0  (nlst) 

W02,rct  =  4p 

(26) 

MH2(n/st) 

WH2,rct  =  2F 

(27) 

...  MH2o(nfst) 

Ww,gen  =  2F 

(28) 

According  to  Ref.  [8],  Wca,out  is  outlet  gas  flow  rate  and  calculated 
as  follows: 

WCa,out  —  kca,out(Pca  ~  Patm ) 

(29) 

where  kca,0ut  is  an  orifice  constant,  Pca  is  gas  pressure  of  cathode 
and  include  the  partial  pressure  of  oxygen  Pq2  ,  the  partial  pressure 
of  nitrogen  PN2  and  the  partial  pressure  of  vapor  Pv,Ca- 

Pca  —  P(J2  +  Pn2  +  Pv,ca 

(30) 

m02Ro2T 

1  °2  W 

vca 

(31) 

mN  PN  T 

Pn2  =  w 

vca 

(32) 

n  mv,caRvT 

1  v,ca  -  w 

vca 

(33) 

Table  1 

Model  parameters  of  PEM  fuel  cell. 


Symbol 

Variable 

Value 

Kh2 

Hydrogen  gas  constant 

4124.3  J  (kg  K)-1 

Rn2 

Nitrogen  gas  constant 

296.8  J  (kg  K)-1 

Ro2 

Oxygen  gas  constant 

259.8 J  (kg  K)"1 

Rv 

Vapor  gas  constant 

461. 5  J  (kgK)-1 

Mh2 

Hydrogen  molar  mass 

2.016  x  10“3  kg  mol”1 

Mq2 

Oxygen  molar  mass 

32  x  10-3  kg  mol-1 

M  n2 

Nitrogen  molar  mass 

28  x  10-3  kg  mol-1 

Ra 

Air  gas  constant 

286.9  J  (kgK)-1 

Patm 

Atmospheric  pressure 

101.325  kPa 

n 

Number  of  cells  in  FC  stack 

381 

Afc 

Fuel  cell  active  area 

280  cm2 

Van 

Anode  volume 

0.005  m3 

Vca 

Cathode  volume 

0.01  m3 

F 

Faraday  number 

96485 

I Patm 

The  humidity  ratio  of  air 

0.5 

Cl 

Constant 

10 

b2 

Constant 

350 

C3 

Constant 

2 

*max 

Maximum  current  density 

2.2 

on  the  75  kW  stacks  used  in  the  FORD  P2000  fuel  cell  prototype 
vehicle  [8,9].  This  stack  consists  of  381  single  fuel  cells,  the  mem¬ 
brane  of  single  fuel  cell  is  Nation  117.  We  linearize  the  FC  system  at 
the  nominal  operating  point  with  stack  current  as  200  A  and  stack 
temperature  as  338  K,  we  can  get  state  space  equation  of  fuel  cell 
system: 

x  =  Aix  +  B\U  ,  . 

y  =  C\X  +  D^u  {  j 

where  x  =  [mH2 ,  mw,an,  m02 ,  mN2 ,  mw?ca]  is  state  variable,  u  =  [  Wiry-, 
WCp]  is  control  input  variable,  y  =  [Vst,  Wca,out]  is  output  variable. 

2.3.  Model  of  PV  system 


The  current-voltage  relationship  at  a  fixed  cell  temperature  and 
solar  radiation  is  expressed  as  follows: 

I  =  IL-  I0[e^v+IRs^a  -  1  ]  -  YjL!!5.  (35) 

Rsh 

where  IL  is  photo  current,  I0  is  the  diode  reverse  saturation  current, 
Rs  is  the  series  resistance,  Rsh  is  the  shunt  resistance,  and  a  is  the 
modified  ideality  factor. 

The  parameters  in  Eq.  (35)  corresponding  to  operation  at 
standard  rating  condition  (SRC)  are  designed:  aref ,  ILref,  Ioref.  To 
determine  the  values  of  these  parameters  [10],  the  three  known 
I-V pairs  at  SRC  are  substituted  into  Eq.  (35). 

For  short  circuit  current:  I  =  Isc,rep  V=  0 

4c, ref  =  4, ref  ~  4, ref  [^ef^ref  -  1 J  -  (36) 

For  open  circuit  voltage:  7  =  0,  V=V0C<ref. 

o  =  4, ref  -  4, ref  [eVoc,ref/aref  -  lj  -  (37) 

At  the  maximum  power  point:  7=7Smp,re/*  V=Vsmp,ref • 

Imp, ref  =  4, ref  ~  4, ref  [e(^/+V,re/'W‘V  _  lj 


Vm 


mp,ref  +  ^mpjef^s 


Rs 


R, 


\ sh 


(38) 


where  Ro2 ,  Rn2  and  Rv  are  gas  constants  of  oxygen,  nitrogen  and 
vapor,  respectively. 

According  to  the  above  nonlinear  model,  there  are  five  state 
variables  x  =  [mn2,  mw^an,  mo2,  m^2,  mW;Ca].  The  parameters  used 
in  this  model  are  given  in  Table  1 .  Most  of  the  parameters  are  based 


Then  the  modified  ideality  factor  a  is  a  linear  function  of  cell 
temperature: 


TcClref 
Tc,  ref 


(39) 
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Table  2 


PV  system’s  model  parameters. 

Symbol 

Variable 

Value 

Rs 

Series  resistance 

0.969  Q 

Rsh 

Shunt  resistance 

199  £2 

$ref 

Total  absorbed  irradiance  at  SRC 

1000  Wm-2 

be, ref 

Short  circuit  current  at  SRC 

4.37  A 

Voc,ref 

Open  circuit  voltage  at  SRC 

42.93  V 

Tc.ref 

Cell  temperature  at  SRC 

298  K 

k 

Boltzmann’s  constant 

1.3805e-23JI<-1 

Q 

Electron  charge 

1.6e-19C 

Table  3 

Model  parameters  of  lithium-ion  battery. 


Symbol 

Variable 

Value 

Ri 

Resistance 

40  x  10“3 

R2 

Resistance 

110 x 10-3 

C 

Capacitance 

4F 

cp 

Specific  heat 

925  J  (kgl<)-1 

m 

Battery  mass 

0.041  kg 

A 

Battery  external  surface  area 

4.3  x  10“3  m2 

Ta 

Ambient  temperature 

23  °C 

where  Tcref  is  the  cell  temperature  for  reference  condition  and  Tc  is 
the  cell  temperature. 

The  diode  reverse  saturation  current  I0  is  calculated  as  follow: 


where  k  is  Boltzmann’s  constant  and  Eg  is  the  material  band  gap. 
Eg  can  be  calculated  as  follows: 

Eg  =  Eg,Tref  [1  -  0.0002677(T  -  Tre/)J  (41 ) 

where  Egjref  =  1.121  eV  for  silicon  cells  and  Egjref  =  1 .6  eV  for  the 
triple  junction  amorphous  cell. 

Photo  current  k  is  calculated  as  follow: 

h  =  —  x  ( ILjef  +  aisc(Tc  -  Tc  ref ))  (42) 

ref 


2.4.  Model  of  lithium-ion  battery 

The  equilibrium  potential  of  the  lithium-ion  battery  depends 
on  the  temperature  and  the  amount  of  active  material  available 
in  the  electrodes,  which  can  be  specified  in  terms  of  state  of  dis¬ 
charge  (SOD).  The  discharge  capacity  of  the  battery  depends  on  the 
discharge  rate  and  the  temperature. 

The  expression  for  the  potential,  the  terminal  voltage  and  the 
SOD  [11]  can  be  given  by 

E[ib(t ),  Tb(t ),  t]  =  v[ib(t ),  Tb{t\  t]  -  R[ntib(t )  (44) 

n 

V[ib(t),  Tb(t),  t]  =  ■  SODkMt),  T„(t),  t]  +  A E(Tb)  (45) 

k= 0 

SOD[ib(t),  Tb(t),  t]  =  Ej  alib(tmTb(t)]  ■  ib(t)dt  (46) 


where  S  is  total  absorbed  irradiance  (Wm-2),  aisc  is  temperature 
coefficient  for  short  circuit  current  (AK-1 ).  Sref  is  total  absorbed 
irradiance  at  SRC  (Wm-2). 

PV  system’s  model  parameters  in  this  paper  are  listed  in  Table  2. 
The  V-I  curve  of  PV’s  single  cell  is  shown  in  Fig.  2.  The  curve 
can  be  divided  into  three  zones:  current  zone,  maximal  power 
zone  and  voltage  zone.  When  0  <  Vpv  <  V\  the  V-I  curve  is  under 
current  zone,  the  state  space  equation  is  Ipv  =  Isc  +  ot\ Vpv.  When 
V\  <  Vpv  <  V2  the  V-I  curve  is  under  maximal  power  zone,  the  state 
space  equation  is  Ipv  =  oi\Vpv  -  a2Vmp  +  /mp.When  V2  <  Vpv  <  Voc 
the  V-I  curve  is  under  voltage  zone,  the  state  space  equation  is 
Ipv  =  a3Vpv  -  a3V0C .  The  parameter  oq,  a2,  a3 ,  Vi,  V2  can  be  calcu¬ 
lated  from  Eqs.  (34)  to  (42). 

State  space  equation  can  be  obtained  by  linearizing  the  PV  sys¬ 
tem: 

x=A2x  +  B2u 

y  =  C2x  +  D2u  K  ; 

where  state  variable  is  Ipv,  control  variable  is  Vpv  and  output  variable 
iS  Jpy 


Fig.  2.  V-I  curve  of  PV’s  single  cell. 


According  to  data  from  a  Sonyl8650  lithium  ion  battery,  a[T(t)] 
and  /3[T(t)]  can  be  calculated  as  follows: 


a[ib{t)]  = 


1.05  if  ib(t)  >  2.8  A 

1.01  + 1.04/1.4  *(ib(t)- 1.4)  if  1.4  <  ib(t)  <  2.8  A 

1  +0.01/0.7  *  (ib(t)  —  0.7)  if  0.7  <  ib(t)  <  1.4  A 

0.96  if  0  <  ib(t)  <  0.3  A 


(47) 


PITM]  = 


1.085  -  0.0015(Tb(t)  +  20)  if  -  20°C  <  Tb(t)  <  -10°C 
1.075  -  0.0007 (Tb(t)  +  10)  if  -  10°C  <  Tb(t )  <  0°C 
1 .063  -  0.00286 Tb{t)  if  -  0 °C  <  Tb(t )  <  22  °C 

1 -0.001(Tft(t)-22)  if  -  22  °C  <  Tb(t)  <  45  °C 


(48) 


A  E[Tb(t)]  = 


-0.46  +  0.012(Tb(f)  +  20)  if  -  20°C  <  Tb(t)  <  -10°C 
-0.38  +  0.025(Tb(t)  +  10)  if  -  10°C  <  Tb(t)  <  0°C 
0.13  +  0.0059Tb(t)  if  0°C  <  Tb(t)  <  22 °C 

0.0104(Tb(t)  —  22)  if  22°C  <  Tb(t)  <  45°C 


(49) 


Combining  the  electrochemical  characteristics  represented 
with  the  equation  for  the  circuit  elements  and  the  circuit  constraint 
equation  yields  the  following  equation  that  relates  the  terminal 
potential  to  the  terminal  current 

W)  =  T  x  [Vb(t)  -  E[ib(t),  Tb(t),  t ]  -  R,ib(t)] 

+  cft  [Vi,(t)  -  E\ib(t),  Tb(t),  t ]  -  R]  i„(t )]  (50) 


The  temperature  change  of  the  battery  is  governed  by  the  ther¬ 
mal  energy  balance  [12]  described  by 

m  ■  cp  ■  =  jb(tf  .  R,  +  T[Vb(t)  -  E[ib(t),  Tb(t),  t]  -  ib(t)Rj  ]2 

-  hcA[Tb(t)  -  Ta]  (51) 


Model  parameters  of  lithium-ion  battery  in  this  paper  are  listed 
in  Table  3. 

State  space  equation  can  be  divided  into  16  sub-models,  where 
ib  is  divided  into  four  zones  according  to  Eq.  (47)  and  Tb  is  also 
divided  into  four  zones  according  to  Eq.  (48).  State  space  equation 
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of  lithium-ion  battery  can  be  obtained  by  linearizing  the  nonlinear 
model  of  lithium-ion  battery, 

x  =  A3x  +  B3u 

y  =  C3x  +  D3u  1  j 


where  state  variable  x  is  [ib,  Tb],  control  variable  u  is  Vb  and  output 
variable  y  is  ib. 


3.  Design  of  model-based  fault  diagnosis 


The  design  of  a  fault  detection  system  requires  obviously  quick 
fault  detection  and  isolation  (FDI)  for  adequate  decision  making. 
Hence,  to  preserve  the  safety  of  operators  and  the  reliability  of  pro¬ 
cesses,  the  presence  of  faults  must  be  taken  into  account  during 
the  system  control  design.  During  the  system  operation,  faults  or 
failures  may  affect  the  sensors,  the  actuators,  or  the  system  com¬ 
ponents.  These  faults  can  occur  as  additive  or  multiplicative  faults 
due  to  a  malfunction  or  equipment  aging. 

The  state  space  representation  of  hybrid  FC  and  photovoltaic  DC 
power  sources  system  is  as  follows: 

y  =  Cx  +  Du  1  j 

where statevariablesxare[mH2,  mw?an,  m02,  mN2,  mw,ca,/Py,  ib ,  Tb], 
output  variables  y  are  [Vst,  Wca,out ,  Ipv,  h]  and  control  variable  u  is 
[Winj,Wcp,Vpv,Vb]. 

After  the  continuous  state  space  model  is  converted  into  discrete 
state  space  model,  the  state  space  representation  of  a  system  that 
may  be  affected  by  actuator  or  sensor  fault  [13]: 

x(k  +  1)  =  Ax{k)  +  Bu{k)  +  Fafa{k)  ,  . 

y{k)  =  Cx{k)  +  Fsfs{k)  ^ 


where  Fa  and  Fb  are  assumed  to  be  known  matrices,  fa  and  fb  are 
the  magnitude  of  the  actuator  and  the  sensor  faults,  respectively. 

In  the  presence  of  sensor  and  actuator  faults,  Eq.  (54)  can  also 
be  represented  by  the  unified  general  formulation: 

f  x(k  +  1 )  =  Ax(k)  +  Bu(k)  +  Fxf{k )  (ss'i 

1  y{k)  =  Cx{k)  +  Fym 


time(second) 


I 


Fig.  3. 

20%. 


timeCsecond) 

Comparison  of  actual  output  and  model  estimation  when  Wca  in  is  decreased 


B  =  U~'B  = 


Bj(k) 

B2(k) 


Then  the  estimation  of  the  actuator  fault  magnitude  is  defined 
as 


/d(fc)  =  WT1 (Xi(fc  +  1 )  -  Auxm  -  AnHk)  ~  (58) 


where  /  =  [fj  fj  ]T  is  a  common  representation  of  sensor’s  and 
actuator’s  faults.  The  fault  vector /in  (55)  can  be  split  into  two 
parts.  The  first  part  contains  the  “d”  faults  to  be  isolated.  In  the 
second  part,  the  other  faults  are  gathered  in  a  vector. 

The  fault  magnitude  estimation  of  the  corrupted  element  is 
extracted  directly  from  the  jth  unknown  input  observer  which 
is  built  to  be  insensitive  to  the  jth  fault  (/*(/<)  =  0).  Based  on  the 
unknown  input  observer,  the  substitution  of  the  state  estimation 
in  the  faulty  description  [13]  leads  to 


Fafdik )  =  x(k  +  1 )  -  Ax(k  +  1 )  -  Bu{k) 


(56) 


In  the  presence  of  an  actuator  fault,  Fd  is  a  matrix  of  full  col¬ 
umn  rank.  Thus,  the  estimation  of  the  fault  magnitude/0(/<)  =fdV<) 
makes  use  of  the  singular-value  decomposition  (SVD)  [14]. 


Let  Fd  =  U 


R 

0 


VT  be  the  SVD  of  Fd.  Thus  R  is  a  diagonal  and 


nonsingular  matrix  and  U  and  V  are  orthonormal  matrices. 
Using  the  SVD  and  substituting  it  into  (56)  and  results  in 


For  sensor  faults,  the  output  equation  is  broken  down  and  can 
be  written  as 


y(k)  =  Cx(k)  +  Esfs(k)  = 


Ci 

C2 


x(k)  + 


Fsi 

FS2 


fs(k) 


(59) 


The  integral  error  vector  z  will  be  affected  by  the  fault  as  well. 
The  integral  error  vector  can  be  described  as  follows: 

z(k+  1)  =  z(/<)  +  Ts(yr(/<)  —  yi(/<)) 

=  z{k)  +  Ts(yr(k)  -  dx(k)  -  Fsi/sW)  1  J 


The  sensor  fault  magnitude  can  be  estimated  as  follows: 

ESXS(1<  +  1 )  =  AsXs(k)  +  BsU(k)  +  Gsyr(k )  (61 ) 


where 


0  0" 
IP  0 
0  Fs 


kk  +  I)=Ak{k)+ACi{k)  + 


VTfd(k) 


where 


A  =  U-'AU  = 


Au(k) 

A2i{k) 


Ai2(/<) 

A22{k) 


(57) 


Xs(k)  = 


"*(/<)" 

z{k) 

lfs(k)  \ 


As  =A 
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Fig.  4.  Time  evolution  of  the  residual  when  Wca  ir,  is  decreased  20%. 


In  order  to  generate  a  set  of  residual,  a  full-order  observer  is 
built  as  follows  [15]: 

f  w(k  +  1)  =  Ew{k)  +  TBu(k)  +  Ky{k)  ,  . 

|  x(/<)  =  w(/<)  +  Hy(k)  1  j 

where  x  is  the  estimated  state  vector  and  w  is  the  state  of  full-order 
observer.  E,  T ,  I<  and  H  are  matrices  to  be  designed  for  achieving 
unknown  input  decoupling  requirements.  The  design  of  E,  T,  I<  and 
H  is  achieved  by  solving  the  following  equations: 


(HC-I)Fd  =  0 

(63a) 

II 

1 

£ 

n 

(63b) 

E  =  A  —  HAC  —  C 

(63c) 

K 2  =  EH 

(63d) 

k  =  ka  +/c2 

(63e) 

If  the  following  conditions  are  fulfilled:  (1  )Rank{CFd)  =  Rank{Fd)1 
(2 )  ( C,  A\ )  is  a  detectable  pair,  where  A\  =  E  +  K\  C,  an  unknown  input 
observer  provides  an  estimation  of  the  state  vector,  used  to  gener¬ 
ate  a  residual  vector: 


(64) 


3P 


model  estimation 

* 

i .   : 

actual  output 

Fig.  5. 
10%. 


timefsecond) 

Comparison  of  actual  output  and  model  estimation  when  Wca  in  is  decreased 


4.  Results  and  discussion 

The  research  results  apply  to  hybrid  power  sources  bus  which 
includes  fuel  cell,  PV  and  lithium-ion  battery.  In  this  hybrid  power 
sources  bus  the  maximal  power  of  fuel  cell  is  75  kW,  the  maxi¬ 
mal  power  of  PV  is  1 .5  kW  and  the  capacity  of  lithium-ion  battery 
is  about  15  Ah.  The  proposed  model-based  fault  diagnosis  method 
is  implemented  in  Matlab/Simulink  environment.  Results  of  three 
of  the  proposed  fault  scenarios  and  f3)  are  presented. /I ,  f2 

and  /3  are  increase  in  the  compressor  motor’s  friction,  humidifier 
choking  and  decrease  in  total  absorbed  irradiance  of  the  PV  system, 
respectively.  According  to  Eq.  (64),  rlt  r2,  r3,  r4  are  the  residuals  of 
[ Vst ,  Wca,out,  Ipv->  ib\- 

As  discussed  in  previous  sections,  the  fault  detection  is  based  on 
checking  the  difference  (residual)  between  the  signals  monitored 
by  a  sensor  and  its  estimation  using  the  detection  model  at  each 
sample  time.  The  fault /I  (increase  in  the  compressor  motor’s  fric¬ 
tion)  is  introduced  into  the  system  at  time  50  s  creating  changes 
in  its  internal  dynamics.  When  the  compressor  motor’s  friction 
increase,  the  mass  flow  rate  of  air  entering  the  cathode  (' Wca  in )  is 
decrease  20%,  the  actual  output  and  model  estimation  of  Vst,  Wca,out 
are  shown  in  Fig.  3(a)  and  (b).  The  residuals  r3,  r4  are  almost  zero 
and  the  residuals  rlt  r2  are  shown  in  Fig.  4(a)  and  (b).  Threshold 
of  Vst  (dash  line  in  Fig.  4(a))  is  about  ±3  V  and  threshold  of  Wca,out 
(dash  line  in  Fig.  4(b)  is  about  ±0.1  kgs-1).  Fig.  4(a)  and  (b)  shows 
that  the  residuals  r\ ,  r2  are  larger  than  their  threshold  after  50  s.  So 
the  residuals  ri,  r2  are  sensitive  to/1  and  the  residuals  r3,  r4  are  not 
sensitive  to  /I .  Then  we  can  detect  the  /I  by  calculating  the  resid¬ 
uals  X\  and  r2.  When  the  compressor  motor’s  friction  increase,  the 
mass  flow  rate  of  air  entering  the  cathode  (Wca  in)  is  decrease  10%, 
the  actual  output  and  model  estimation  of  Vstl  Wca,out  are  shown  in 
Fig.  5(a)  and  (b).  The  residuals  r3,  r4  are  almost  zero  and  the  resid¬ 
uals  ri,  r2  are  shown  in  Fig.  6(a)  and  (b).  Fig.  6(a)  and  (b)  shows 
that  the  residuals  r\ ,  r2  are  smaller  than  their  threshold  after  50  s 


r(k)=y(k)-Cx(k) 


((kgs'1)  cr  F„.(V)  “  r2(kgs-')  cr  r,(V) 
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Time  evolution  of  the  residual  when  Wcatin  is  decreased  10%. 
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Fig.  8.  Time  evolution  of  the  residual  corresponding  to  fault  fl. 
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and  the  fault /I  is  not  detected.  So  the  result  of  fault  detection  is 
sensitive  to  the  magnitude  of  fault. 

The  fault  fl  is  the  choking  of  humidifier.  The  fault  f2  is  intro¬ 
duced  into  the  system  at  time  50  s  creating  changes  in  its  internal 
dynamics.  When  the  fault  f2  is  applied,  the  actual  output  and  model 
estimation  of  Vst,  Wca,out  are  shown  in  Fig.  7(a)  and  (b).  The  residuals 
r3,  r4  are  almost  zero  and  the  residuals  rlt  r2  are  shown  in  Fig.  8(a) 
and  (b).  Fig.  8(a)  and  (b)  shows  that  the  residual  r\  is  larger  than 
the  threshold  and  r2  is  less  than  the  threshold  after  50  s.  This  shows 
that  fault  f3  appear.  So  the  residual  r\  is  sensitive  to  f2  and  the  resid¬ 
uals  r2,  r3,  r4  are  not  sensitive  to  fl.  Then  we  can  detect  the  fl  by 
calculating  the  residuals  r2. 

The  fault  f3  is  decrease  in  total  absorbed  irradiance  of  the  PV  sys¬ 
tem.  The  fault  f3  is  introduced  into  the  system  at  time  50  s  creating 
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Fig.  7.  Comparison  of  actual  output  and  model  estimation  when  fl. 
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Fig.  9.  Comparison  of  actual  output  and  model  estimation  corresponding  to  fault  f3. 
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Fig.  10.  Time  evolution  of  the  residual  corresponding  to  fault  f3. 

changes  in  its  internal  dynamics.  When  the  fault  /3  is  applied,  the 
actual  output  and  model  estimation  of  Ipv  are  shown  in  Fig.  9.  The 
residuals  rlt  r2,  r4  are  almost  zero  and  the  residual  r3  is  shown  in 
Fig.  1 0.  Threshold  of  lpv  (dash  line  in  Fig.  1 0)  is  about  ±1 .5  A.  Fig.  1 0 
shows  that  the  residual  r3  is  larger  than  the  threshold  after  50  s. 
This  shows  that  fault  f3  appear.  So  the  residual  r3  is  sensitive  to  f3 
and  the  residuals  r\ ,  r2,  r4  are  not  sensitive  to  f3.  Then  we  can  detect 
the  J3  by  calculating  the  residuals  r3. 

From  the  simulation  results  we  can  see  that  we  can  detect  fault 
fl,J2  and  f3  by  calculating  the  corresponding  residuals  rlt  r2  and  r3. 
This  shows  that  the  proposed  model-based  fault  detection  method 
can  detect  faults  of  hybrid  DC  power  sources  on-line. 

5.  Conclusions 

In  this  paper,  a  new  model-based  fault  detection  methodology 
of  hybrid  DC  power  sources  has  been  presented  and  tested.  An 
advantage  of  this  new  methodology  is  that  it  can  find  the  fault 
on  line,  improve  the  generation  time  and  avoid  permanent  dam¬ 
age  to  the  equipment.  To  prove  this  methodology,  a  hybrid  DC 
power  sources  Matlab  simulator  based  on  the  state  space  model 
is  built.  The  simulator  is  modified  to  include  a  set  of  possible  fault 


scenarios  proposed  in  this  work.  This  modified  simulator  allows 
imposing  a  determined  fault  scenario,  within  the  considered  set 
of  faults  in  the  hybrid  DC  power  sources  and  analysis  its  behavior. 
All  the  simulated  faults  have  been  tested  with  the  new  fault  detec¬ 
tion  methodology.  The  results  show  that  the  proposed  model-based 
fault  detection  method  can  detect  faults  of  hybrid  DC  power  sources 
effectively. 
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